library("rstudioapi")     
##Set working directory
setwd(dirname(getActiveDocumentContext()$path))
getwd()
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(xlsx)
library(estimatr)
library(tidyverse)
library(geosphere)
library(DescTools)
require(dplyr)

load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

working_final = final_finalV2

# Cobalt package and related models: 
require(cobalt)

# rename variables:  
working_final$HHSize = working_final$Q141
working_final$NonFood = working_final$Q145
working_final$Employment = working_final$Q143
working_final$Employment = as.factor(working_final$Employment)
working_final$Children = working_final$Q142
working_final$Finance = working_final$Q146
working_final$Gender = as.factor(working_final$Gender)

file1 = working_final[,c("individual_treatment", 
                         "Age","Gender","HHSize",
                         "Children","Employment","Income", 
                         "NonFood","Finance", "Attended")]

dim(file1)
dim(na.omit(file1))
file1 = na.omit(file1)

# Bal Tab no. 1: 
# Balance summary across all treatment pairs
bal.tab(individual_treatment ~ 
          Age + Gender + HHSize +
          Children + Employment + Income + 
          NonFood + Finance + Attended, 
          data = file1, 
          weights = "att.weights", 
          binary = "std", continuous = "std") 


# Balance tables with mean differences, variance ratios, and
# statistics for the unadjusted sample
bal.tab(individual_treatment ~ Age + Gender + HHSize +  Children + Employment + Income + 
        NonFood + Finance + Attended,
        data = file1,
        weights = "att.weights",
        disp = c("means", "sds"), 
        un = TRUE,
        stats = c("mean.diffs", "variance.ratios"))



#Generating propensity score weights for the ATT
library(WeightIt)
W.out <- WeightIt::weightit(individual_treatment ~ Age + Gender + HHSize +  Children + Employment + Income + 
                            NonFood + Finance + Attended,
                            data = file1,
                            method = "ps")
bal.tab(W.out)


#Using WeightIt to generate weights with multinomial
#logistic regression
#Max.Diff.Un Max.Diff.Adj
W.out.mn <- WeightIt::weightit(individual_treatment ~ Age + Gender + HHSize +  Children + Employment + Income + 
                               NonFood + Finance + Attended,
                               data = file1,
                               method = "ps",
                               use.mlogit = FALSE)
#Balance summary across treatment pairs
bal.tab(W.out.mn, un = TRUE)

#Assessing balance for each pair of treatments
# with means:
bal.tab(W.out.mn, un = TRUE,
        disp.means = TRUE,
        which.treat = .all)


# Most important graph: 
#Summarizing balance in a Love plot
love.plot(W.out.mn, thresholds = c(m = .1), binary = "std",
          which.treat = .all, abs = FALSE)

pdf("ExtData_Figure1.pdf", width = 12.5, height =9.5)
love.plot(W.out.mn, thresholds = c(m = .1), binary = "std",
          which.treat = .all, abs = FALSE)
dev.off()



